% The actual Travel model and modifications
%
% Input:  PreferenceFactor (used for minimization),cfg,exvars,depvars,indepvars,r,t
% Output: Residuals (depeding on calibration type; only during calibrations)
%         depvars structure
%
function varargout = ModelTravelLinearTW(cfg,exvars,depvars,indepvars,r,time,varargin)
%% Destinguish between calibration or validation and specific type of run
    for n = 2:2:length(varargin)
         eval( [varargin{n-1} '(r,time(1),:)=varargin{n};']);   
    end
    
    modes                             = indepvars.UseModes(r,1,:); % get rid of modes not used by Schafer (e.g. High speed train)
    if cfg.ConvergencePFandSW && sum(r==[1 2 3])==0 %&& isempty(varargin);
        
        depvars.PreferenceFactor(r,:,:)   = repmat(indepvars.PreferenceFactorStart(r,:,:), [1 length(cfg.Time) 1]) ...
        + repmat(cfg.Time,size(indepvars.PreferenceFactorStart(r,:,:))) .* ...
        repmat(indepvars.PreferenceFactorSlope(r,:,:),[1 length(cfg.Time) 1]); % This definition allows for a slope in the PF
         
        depvars.PreferenceFactor(r,:,3)   =  exvars.Income(r,:) ./((sum(exvars.Income(1:3,:).*exvars.Population(1:3,:),1))./sum(exvars.Population(1:3,:),1)).*indepvars.PFSaturation(r);    
    else

        depvars.PreferenceFactor(r,:,:)   = repmat(indepvars.PreferenceFactorStart(r,:,:), [1 length(cfg.Time) 1]) ...
        + repmat(cfg.Time,size(indepvars.PreferenceFactorStart(r,:,:))) .* ...
        repmat(indepvars.PreferenceFactorSlope(r,:,:),[1 length(cfg.Time) 1]); % This definition allows for a slope in the PF
    end
         
    indepvars.TravelPrice(r,:,:) = exvars.TravelPrice(r,:,:);
    if indepvars.ComfElasticity(r) ~= 0.5
    exvars.EnergyCost(exvars.EnergyCost<0) =0;
    indepvars.TravelPrice(r,35:130,:) = ( exvars.NonEnergyCost(r,35:130,:) ...
        .*((repmat(exvars.Income(r,35:130,:), [1 1 size(depvars.TravelDemand,3)])...
        ./repmat(exvars.Income(r,35,:), [1 96 size(depvars.TravelDemand,3)]))...
        .^-0.5...
        .*(repmat(exvars.Income(r,35:130,:), [1 1 size(depvars.TravelDemand,3)])...
        ./repmat(exvars.Income(r,35,:), [1 96 size(depvars.TravelDemand,3)]))...
        .^indepvars.ComfElasticity(r,1,1)) ...
        + exvars.EnergyCost(r,35:130,:)) ...
        ./ exvars.LoadFactor(r,35:130,:);
    end
    
for t = time
    for n = 1:indepvars.TimeWeightRuns % Aditional iterations make TimeWeight function better
%% Girod et. al TIMER model
    if t==1 && n==1% checks if timeweight is a calibration variable
    depvars.TimeWeight(r,1)  = indepvars.TimeWeight(r,1);
    elseif n==1
        depvars.TimeWeight(r,t)  =...
        depvars.TimeWeight(r,t-1) .* ...
        sum( depvars.TravelDemand(r,t-1,modes )./exvars.Speed(r,t-1,modes),3) ...
        ./ indepvars.TTB(r,t);
    else
        depvars.TimeWeight(r,t)  =...
        depvars.TimeWeight(r,t) .* ...
        sum( depvars.TravelDemand(r,t,modes )./exvars.Speed(r,t,modes),3) ...
        ./ indepvars.TTB(r,t);
    end
    
    depvars.TimeWeight(r,t) = max([depvars.TimeWeight(r,t), 0.1]); % Ensures timeweight does not become too small

    depvars.RelativeModeCostTravel(r,t,modes ) =...
    indepvars.TravelPrice(r,t,modes ) ./ (mean(indepvars.TravelPrice(r,t,modes )) .* depvars.PreferenceFactor(r,t,modes))...
    + 1./exvars.Speed(r,t,modes) ./ mean(1./exvars.Speed(r,t,modes)) .* depvars.TimeWeight(r,t);

    depvars.ModeSplitOptimum(r,t,modes ) =...
    exp(-indepvars.Lambda(r) .* depvars.RelativeModeCostTravel(r,t,modes )) ...
    / sum(exp(-indepvars.Lambda(r) .* depvars.RelativeModeCostTravel(r,t,modes )));
        
    if  t==1 % checks if timeweight is a calibration variable
        depvars.ModeSplit(r,t,modes ) = depvars.ModeSplitOptimum(r,t,modes );
    else
        depvars.ModeSplit(r,t,modes )= ...
        depvars.ModeSplit(r,t-1,modes ) + indepvars.Inertia(r,1,modes ) .*( depvars.ModeSplitOptimum(r,t,modes ) - depvars.ModeSplit(r,t-1,modes ) ) ;        
    end
    
    depvars.TMB(r,t)                = (exvars.TMBCorrection(r,1) + 0.035 + 0.085 * sum(depvars.ModeSplit(r,t,cfg.ModesFasterThanLDV),3)) + indepvars.TMBAdded(r); 

    depvars.TravelDemand(r,t,modes) = ...
       depvars.ModeSplit(r,t,modes ) .* depvars.TMB(r,t) .* exvars.Income (r,t) ./ sum(depvars.ModeSplit(r,t,modes ) .* indepvars.TravelPrice(r,t,modes),3);

%% Determine the errors and objective functions

    if t<36 % ensures function is only called with historic data
        depvars.ModeWeight(r,t,modes) = sqrt(exvars.TravelDemand(r,t,modes))./sum(sqrt(exvars.TravelDemand(r,t,modes)),3);% Ensures not all modes are equal, but that car is not the only important mode
        depvars.ObjectiveFunction(r,t) = sqrt(sum((depvars.TravelDemand(r,t,modes) - exvars.TravelDemand(r,t,modes)).^2 ...
        ./ exvars.TravelDemand(r,t,modes).^2 .* depvars.ModeWeight(r,t,modes).^2,3)); % used for optimization
%           depvars.ObjectiveFunction(r,t) = sqrt((( sum(depvars.TravelDemand(r,t,modes) - exvars.TravelDemand(r,t,modes) ,3).^2) ...
%          ./sum(exvars.TravelDemand(r,t,modes),3).^2));
        depvars.Resid(r,t)             =  depvars.ObjectiveFunction(r,t);
    end
    
    end
end
%% Write output
    varargout{1} = mean(depvars.ObjectiveFunction(r,time));
if nargout==2
    varargout{2} = depvars;
end
    
%% Visualizes minimization process
if cfg.LiveMinimationVisual == 1 && ismember(time(1),cfg.LiveMinimizeVisuals{2}); % bit clumsy to make the code faster
     GeneralPlotting('LiveMinimizeVisuals','Dont save',cfg,exvars,depvars,indepvars,r,time);
elseif cfg.PlotMinimize == 1
    GeneralPlotting('SaveToGlobal',' ',cfg,exvars,depvars,indepvars,r,time);
end
    
end